ConvertGeodeticToTransverseMercator Subroutine

SUBROUTINE ConvertGeodeticToTransverseMercator &
!
(lon, lat, k, centM, lat0, a, e, eb, falseN, falseE, x, y)

USE Units, ONLY: &
!Imported parametes:
pi

USE StringManipulation, ONLY: &
!Imported routines:
ToString

IMPLICIT NONE

! Arguments with intent(in):
REAL (KIND = float), INTENT (INOUT) :: lon !!geodetic longitude [radians]
REAL (KIND = float), INTENT (IN) :: lat !!geodetic latitude [radians]
REAL (KIND = float), INTENT (IN) :: k !!scale factor
REAL (KIND = float), INTENT (IN) :: centM !!central meridian [radians]
REAL (KIND = float), INTENT (IN) :: lat0 !!latitude of origin [radians]
REAL (KIND = float), INTENT (IN) :: a !! semimajor axis [m]
REAL (KIND = float), INTENT (IN) :: e !! eccentricity
REAL (KIND = float), INTENT (IN) :: eb !! second eccentricity
REAL (KIND = float), INTENT (IN) :: falseN !!false northing
REAL (KIND = float), INTENT (IN) :: falseE !!false easting



!Arguments with intent(out):
REAL (KIND = float), INTENT (OUT) :: x !!easting coordinate [m]
REAL (KIND = float), INTENT (OUT) :: y !!northing coordinate [m]

!Local variables:
REAL (KIND = float) :: c !Cosine of latitude
REAL (KIND = float) :: c2
REAL (KIND = float) :: c3
REAL (KIND = float) :: c5
REAL (KIND = float) :: c7
REAL (KIND = float) :: es !!eccentricity squared
REAL (KIND = float) :: ebs !!second eccentricity squared
REAL (KIND = float) :: dlam  !Delta longitude - Difference in Longitude
REAL (KIND = float) :: eta !constant - TranMerc_ebs *c *c
REAL (KIND = float) :: eta2
REAL (KIND = float) :: eta3
REAL (KIND = float) :: eta4
REAL (KIND = float) :: s !Sine of latitude
REAL (KIND = float) :: sn !Radius of curvature in the prime vertical
REAL (KIND = float) :: t !Tangent of latitude
REAL (KIND = float) :: tan2
REAL (KIND = float) :: tan3
REAL (KIND = float) :: tan4
REAL (KIND = float) :: tan5
REAL (KIND = float) :: tan6
REAL (KIND = float) :: t1 !Term in coordinate conversion formula
REAL (KIND = float) :: t2 !Term in coordinate conversion formula
REAL (KIND = float) :: t3 !Term in coordinate conversion formula
REAL (KIND = float) :: t4 !Term in coordinate conversion formula
REAL (KIND = float) :: t5 !Term in coordinate conversion formula
REAL (KIND = float) :: t6 !Term in coordinate conversion formula
REAL (KIND = float) :: t7 !Term in coordinate conversion formula
REAL (KIND = float) :: t8 !Term in coordinate conversion formula
REAL (KIND = float) :: t9 !Term in coordinate conversion formula
REAL (KIND = float) :: tmd !True Meridional distance
REAL (KIND = float) :: tmdo !True Meridional distance for latitude of origin
REAL (KIND = float) :: temp_Origin
REAL (KIND = float) :: temp_Long

! Local parameters:
REAL (KIND = float), PARAMETER :: MAX_LAT = 89.99 * degToRad    ! 89.99 degrees in radians
REAL (KIND = float), PARAMETER :: MAX_DELTA_LONG =  90. * degToRad    ! 90. degrees in radians

!------------end of declaration------------------------------------------------

IF ( lat < -MAX_LAT .OR. lat > MAX_LAT ) THEN
  CALL Catch ('error', 'GeoLib',   &
    'Converting Geodetic to Transverse Mercator: &
    latitude out of range' ,  &
    code = consistencyError, argument = ToString(lat*radToDeg)//' deg' )
END IF
  
IF ( lon > pi) THEN
   lon = lon - 2*pi
END IF

 IF ( lon < centM - MAX_DELTA_LONG .OR. lon > centM + MAX_DELTA_LONG ) THEN
    IF ( lon < 0. ) THEN
      temp_Long = lon + 2 * pi
    ELSE
      temp_Long = lon
    END IF
    
    IF ( centM < 0. ) THEN
      temp_Origin = centM + 2 * pi
    ELSE
      temp_Origin = centM
    END IF
    
    IF ( temp_Long < temp_Origin - MAX_DELTA_LONG .OR. &
         temp_Long > temp_Origin + MAX_DELTA_LONG   ) THEN
        CALL Catch ('error', 'GeoLib',   &
         'Converting Geodetic to Transverse Mercator: &
         longitude out of range' ,  &
         code = consistencyError, argument = ToString(lon*radToDeg)//' deg' ) 
    END IF
    
 END IF
 
 !Delta longitude
 dlam = lon - centM
 
 ! Check for distortion and send warning.
!Distortion will result if Longitude is more than 9 degrees 
!from the Central Meridian at the equator
!and decreases to 0 degrees at the poles
!As you move towards the poles, distortion will become more significant 
IF ( ABS(dlam) > (9.0 * degToRad) * COS(lat) ) THEN
   CALL Catch ('warning', 'GeoLib',   &
     'Converting Geodetic to Transverse Mercator: &
     possible distortion in longitude computation ' ,  &
      argument = ToString(lon*radToDeg)//' deg' ) 
END IF

IF ( dlam > pi ) THEN
  dlam = dlam - (2 * pi)
END IF

IF ( dlam < -pi ) THEN
  dlam = dlam + (2 * pi)
END IF

IF ( ABS(dlam) < 2.e-10 ) THEN
  dlam = 0.0
END IF

!Eccentricities squared
es = e**2.
ebs = eb**2.

! Sine Cosine terms
s = SIN (lat)
c = COS (lat)
c2 = c * c
c3 = c2 * c
c5 = c3 * c2
c7 = c5 * c2

! Tangent Value
t = TAN (lat)
tan2 = t * t
tan3 = tan2 * t
tan4 = tan3 * t
tan5 = tan4 * t
tan6 = tan5 * t
eta = ebs * c2
eta2 = eta * eta
eta3 = eta2 * eta
eta4 = eta3 * eta

! Radius of curvature in the prime vertical
sn = Nu (lat, a, e)

!True Meridianal Distances 
tmd = MeridionalDistance (lat, a, e)  

! True Meridional Distance for latitude of origin 
tmdo = MeridionalDistance (lat0, a, e)

!northing 
t1 = (tmd - tmdo) * k
t2 = sn * s * c * k/ 2.e0
t3 = sn * s * c3 * k * (5.e0 - tan2 + 9.e0 * eta + 4.e0 * eta2) /24.e0
t4 = sn * s * c5 * k * (61.e0 - 58.e0 * tan2                     &
     + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2 &
     + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4        &
     -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0

t5 = sn * s * c7 * k * (1385.e0 - 3111.e0 * &
     tan2 + 543.e0 * tan4 - tan6) / 40320.e0

y = falseN + t1 + dlam**2. * t2 + dlam**4. * t3 + dlam**6. * t4 + dlam**8. * t5

!Easting 
t6 = sn * c * k
t7 = sn * c3 * k * (1.e0 - tan2 + eta ) / 6.e0
t8 = sn * c5 * k * (5.e0 - 18.e0 * tan2 + tan4                       &
     + 14.e0 * eta - 58.e0 * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3 &
     - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 ) / 120.e0
t9 = sn * c7 * k * ( 61.e0 - 479.e0 * tan2 &
     + 179.e0 * tan4 - tan6 ) / 5040.e0

x = falseE + dlam * t6 + dlam**3. * t7 + dlam**5 * t8 + dlam**7. * t9 

END SUBROUTINE ConvertGeodeticToTransverseMercator